Noise driven unlimited population growth 
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Demographic noise causes unlimited population growth in a broad class of models which, without 
noise, would predict a stable finite population. We study this effect on the example of a stochastic 
birth-death model which includes immigration, binary reproduction and death. The unlimited 
population growth proceeds as an exponentially slow decay of a metastable probability distribution 
(MPD) of the population. We develop a systematic WKB theory, complemented by the van Kampen 
system size expansion, for the MPD and for the decay time. Important signatures of the MPD is 
a power-law tail (such that all the distribution moments, except the zeroth one, diverge) and the 
presence in the solution of two different WKB modes. 

PACS numbers: 87.23.Cc, 02.50.Ga 



Since the celebrated essay of Malthus [l| quantitative 
modeling of population dynamics has attracted much in- 
terest. To a large extent, this interest is powered by the 
danger of a Malthusian catastrophe, when too a rapid 
population growth causes a fatal lack of resources. Here 
we focus on a variant of Malthusian catastrophe by con- 
sidering not too a large population that undergoes bi- 
nary reproduction, immigration and death. Although 
macroscopically stable, this population can be pushed 
to the Malthusian limit (a critical population size that 
sparks a Malthusian catastrophe) by rare large fluctua- 
tions. We show that unlimited population growth pro- 
ceeds as a slow decay of a metastable probability distribu- 
tion (MPD) of the population. We determine the MPD 
and the decay time analytically by developing a system- 
atic WKB theory for the master equation and combining 
it with the van Kampen system size expansion. We show 
that the MPD is described by two different WKB modes 
which are strongly coupled in a narrow region around the 
unstable fixed point of the deterministic rate equation of 
the model. At large population sizes the MPD exhibits 
a power-law tail so that all the distribution moments, 
except the zeroth one, diverge. 

Deterministic rate equation. At the deterministic level 
of modeling a Malthusian catastrophe does not occur if 
the gain and loss processes balance each other so that 
the resulting steady-state population size is stable with 
respect to small perturbations. Real populations, how- 
ever, behave stochastically, rather then deterministically 
The stochasticity may cause an unlimited popula- 
tion growth in a broad class of models which determin- 
istic counterparts predict a stable finite population size. 
We will investigate this remarkable phenomenon on the 
example of a birth-death model [3j which accounts for 

binary reproduction 2A — > 3A, immigration A, and 
death A — > 0. The rate equation for this model is 

h = a - fin+ (X/2)n 2 , (1) 

where n(t) 3> 1 is the average population size. For a 
relatively low death rate, /i 2 < 2aX, Eq. (JTJ does not 



have fixed points, and the population size blows up in 
finite time for any fi{t — 0). For ^ 2 > 2a X Eq. (Q} has 
two fixed points: n\ = 0(1 — S) and n% = Q(l + S), where 
= fi/X ^> 1 and S 2 = 1 — 2aX/[i 2 . When starting 
from any n(t = 0) < ri2, the population size flows to 
the attracting fixed point n = m with a characteristic 
relaxation time t t = 1/(jj,S), and stays there forever. 

Master equation, absorbing state and decay of 
metastable state. The demographic noise, ignored by the 
rate equation ([1]), is accounted for by the master equa- 
tion, see e.g. Ref. which governs the evolution of 
probability P n (t) to have n individuals at time t: 



P-n — A n _i-P n _i — (A n + it n )Pn + Mn+i-fn+i ■ 



(2) 



where A„ = (A/2) n(n— 1)+<t, and \x n = \xn. One striking 
property of Eq. {2} concerns its steady state. Summing 
up the first n equations in Eq. Ip]). we obtain 



1=0 



!^ PnHn+1)Pn+1 -2fp n , (3) 



where 7 = 1 — S 2 = 
the death rate: /i£ 



2crA//i 2 , and the time is rescaled by 
-* t. Putting d/dt = 0, we obtain 



P 



n+l 



n(n - 1) +7fi 2 
2fi(n+ 1) 



P, 



(4) 



Clearly, lim P n — 00 unless P n — 0, n — 0,1,2,.... 

n — >oo 

Therefore, the ultimate state of the stochastic process 
corresponds to an empty system, in a stark contrast to 
the prediction of the rate equation ([1]). At the same 
time, there is a separate absorbing state at infinity which 
"collects" the individuals and ultimately becomes fully 
populated, as its probability 7' 00 (i — * 00) = 1. Here 
is an overview of how the unlimited population growth 
occurs. At t > r r , the MPD sets in, peaked at the at- 
tracting fixed point n = n\. We will assume (and check a 
posteriori) that the width of the MPD here is much less 
than the distance n% — n\ between the two fixed points 
of the rate equation. In this regime a large fluctuation 
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is needed to bring the population beyond the unstable 
fixed point n = n 2 of the rate equation, from where it 
rapidly escapes to infinity. As large fluctuations occur 
with certainty, a full transfer of the population to infinity 
is certain. At t » r r (in the physical units), P n (t) de- 
cays as P n (t) ~ Cn n exp(— t/r), whereas Voo(t) grows as 
Vooit) — 1— Ccxp(— i/r). Here 7r„ is the quasi-stationary 
probability distribution (QSD; it is normalized to unity), 
and C is a constant depending on the initial condition: 
for "macroscopic" initial conditions C ~ 1. As in other 
instances of weak-noise-driven escape from a metastable 
state [4]. the decay time r turns out to be exponentially 
long compared to the relaxation time r r . 

The QSD n n and the decay time r are determined by 
the eigenvalue problem 

- Eir n = — [(n — l)(n - 2)n n ^ 1 - n(n - l)n n ] 

+ [(n + l)7T„ + i - mr n ] + — (7T n _i - 7T„) (5) 

where E = {^t)^ 1 > is the rescaled eigenvalue, and 
we are interested in the eigenmode with the smallest E. 
We will exploit the large parameter O 3> 1 and solve the 
eigenvalue problem analytically by combining a system- 
atic WKB expansion [5| with the van Kampen system 
size expansion Q ■ 

The WKB analysis of the master equation that 
we develop here, extends the existing approaches [1, 0, H, 
7], as it accounts for two different WKB modes, and for 
mode-coupling effects, see below. The WKB ansatz is 



a(n) e 



-S(n) 



(6) 



where, for n > 1, we can treat the action 5(n) and am- 
plitude a(n) as continuous functions of n. As can be 
checked a posteriori, the ordering of terms is the fol- 
lowing: 5(n) = 0(0), a{n) = 0(1), S'(n) = 0(1), 
a'(n) = 0(1/0), 5"(n) = 0(1/0), a"{n) = 0(1/O 2 ), 
etc. Here and in the following the primes stand for n- 
derivatives. Therefore, we can approximate 



7Tn±l 



S" a 1 
1 ± — 

2 a 



(7) 



Now we substitute Eqs. © and in Eq. ©. As E 
turns out to be exponentially small in 1/0, we must put 
E = in all WKB orders. In the leading order we obtain 
the eikonal equation H(n,p) = which describes the 
trajectories of the time-independent Hamiltonian 



IhH.p) - (,"- 



7 

~2 



■V ■ (8) 



Here n is the coordinate, and p = S' is the momentum 
0. The zero-energy lines, 

p = Ps =0 and p =p / = -ln(— + — j , (9) 



describe the slow and the fast (as functions of n) WKB 
modes, respectively. For the slow mode the action 5 = 0. 
One can check that the corresponding Hamilton equation 
for n coincides with the rate equation ([T]). The fast mode 
corresponds to the instanton: a heteroclinic orbit of the 
Hamiltonian ([5]) which exits the saddle point (ni,0) and 
enters the saddle point (ri2,0) of the phase plane (n,p), 
see Fig.J^l In analogy with other problems of noise driven 
escape 0] , this instanton describes the most probable 
escape path: in this case to infinity. The action S(n) = 
J n Pf(n')dn' along the instanton is 

n in ^0 \ 

S(n) = n 20^7 arctan — - n In ^— + — j , 

(10) 

where the integration constant can be put to zero. Note 
that the saddle points (m, 0) and (712, 0) of the Hamilto- 
nian {SJ are mode- crossing points, as pf = p s = there. 
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FIG. 1: (color online). The zero-energy trajectories (|9]) on 
the phase plane (n,p) for 7 = 1/2. The thick line indicates 
the WKB modes contributing to the QSD. 



In the subleading order of the WKB expansion we ob- 
tain a first-order equation for the amplitude a[n): 



[a 2 (n 2 eP - 2nfle- p + 7OV)]' 
= 2 (Qe-P - 2neP + n) a 2 , 

For the fast mode we find, after some algebra, 



a f (n) = 



y/n~{n 2 +7O 2 ) ' 



where Af = const 



(11) 



(12) 



Therefore, the fast WKB mode is well behaved. For the 
slow mode Eq. (fTT|) yields 



a s (n) 



(n - ni)(n - n 2 ) 



where A s = const . (13) 



The slow- mode solution diverges at each of the two mode- 
crossing points implying breakdown of the WKB approx- 
imation there. To understand the mechanism of break- 
down, we notice that it occurs in the regions of small 
p = 5', that is a slow variation of 5(n) and, therefore, 
of 7r„. Here we can use the (stationary) Fokker-Planck 
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equation which follows from the van Kampen system size 
expansion applied to Eq. © with E = (ioj: 

[(n -mXna-n) *■„]'+ [(n + ft) 2 - <5 2 ft 2 ] 7r n | =0. 

(14) 

The first and second terms describes drift and diffusion, 
respectively. The mechanism of breakdown of WKB be- 
comes clear once we observe that the slow mode solution 
tt„ = a s (n), as described by Eq. (fl3|) . solves Eq. (fl4|) 
with the diffusion term neglected. As we will see shortly, 
an account of the small diffusion term regularizes the sin- 
gularity. This regularization is needed only in a narrow 
boundary layer around the mode-crossing point n = n<i. 
Indeed, at 1 <C n < n 2 the slow-mode solution is merely 
an exponentially small correction to the fast-mode so- 
lution, and so it should be discarded there. The situ- 



ation is different at 



> 



n 2 . 



Here the fast-mode solu- 



tion af(n)e~ st - n ^ should be discarded, as it diverges as 
n — > oo, whereas the slow mode yields the correct solu- 
tion. The slow and fast modes are strongly coupled in 
the boundary layer around n = n 2 , and this coupling is 
described by the boundary-layer solution. 

The boundary layer and asymptotic matching. Let us 
consider the stationary Fokker-Planck equation (TT4|) in 
the vicinity of the mode-crossing point n — n 2 : \n— n 2 \ <C 
n 2 . Here we can put n — n\ ~ n 2 — n± in the drift 
term, and n ~ n 2 in the diffusion term. Integrating the 
equation once, we obtain 



dn(x)/dx — 2xtt(x) = —C\ 



1/2 

where x = (n — n 2 )/l 2 , l 2 = [20(1/5 + 1)] y is the char- 
acteristic width of the boundary layer, and G\ > is 
a constant. The general solution of Eq. (fl~5|) is ir(x) — 
C\(j){x) +C 2 e x , where <f)(x) — e x e~^ d£, and C 2 is 
another constant which, as we will see shortly, must be 
put to zero. The function <p(x) has the following asymp- 
totes at \x\ ^> 1: 



_ f (2X)- 1 + O (x~ 3 ) 
) — 1 . rz^x 2 , n /ui-i 



1, 



+ 0(\x\~ 1 ) , x < 0, -x> 1 



(16) 



We start the matching procedure from the region of 
n > n 2 , where the solution can only include the slow 
mode (I13p with a yet unknown normalization constant 
A s . Consider the region < n — n 2 -C n 2 . In the leading 
order, Eq. (fT3"|) yields 7r„ = a s (n) ~ A s /[2Q5(n — n 2 )\. 
Matching this asymptote with the boundary layer solu- 
tion tt(x) in their joint region of validity l 2 -C n — n 2 <C 
n 2 , we obtain C\ — A s {£ll?,$)~ x and C 2 = 0. Having 
found C\, we have determined, up to A s , the boundary 
layer solution w(x). Now we match this solution with 
the fast-mode WKB solution a/(n)e~ s (") in their joint 
region of validity l 2 <^ n 2 — n <ti n 2 . To this end we ex- 
pand S(n) from Eq. (|10|) around n — n 2 up to (n — n 2 ) 2 



and evaluate a/(n), given by Eq. (|12[) , at n = n 2 . The 
matching yields 



A f = A^tt/S) 1 / 2 ^! + S)e s ^ 



(17) 



and determines, up to A s , the complete WKB solution 
at 1 -C n < n 2 . 

The WKB approximation breaks down at n — 0(1). 
To find the QSD in this region we return to Eq. ([5]) and 
notice that, at n <C jQ, 7r ra grows rapidly with n, so 
that 7r n _i -C 7r n . The leading-order terms here are the 
following: (n + l)7r n +i — ('-fD,/2)ir n ~ 0. That is, immi- 
gration and death balance each other and dominate over 
the reproduction. The resulting recursion relation yields 
a Poisson distribution: 



7R) 

n! 



2«y 

2 



(18) 



To determine the unknown constant ttq we can match the 
asymptote (TT51) with the asymptote of the WKB solution 
a/(n)e _s< ' n ' at 1 < n < 7ft. To this end we expand the 
action S(n) at n <C 7ft: S(n) ~ — n + nln [2n/(7f2)]. On 
the other hand, at n > 1 we can use Stirling's formula 
nl ~ \j2im (n/e)" in Eq. I|18p. The matching yields 



tto 



2ttA s e s (™ 2 ) 
0^/2(1 -5) 



(19) 



By now we have found, up to the normalization constant 
A s , the QSD for all n. The normalization, in the lead- 
ing order, is determined by the region of |n — n±\ <C ri\, 



(15) where the fast-mode solution 7r n = af(n) 



,-S(«) 



exp 



-(n - ni) 2 



I 2 
'1 



proximately gaussian: 

^A s (l + S)e ilAs 

^ ~ V2 3 / 2 (5 1 /2(l_ 5)3/2 

1 /2 

Here we have denoted h = [2Q(l/5 - 1)] 1 and 
As(tf) = [S(n 2 ) - S(nx)]/tt = 26 

—2\J 1 — S 2 ( arctan 



is ap- 



(20) 



1-5 



arctan 



1-5 



Normalizing the gaussian (|20[) to unity, we obtain 

_ ns(i-5) QAs 

s ~ n(l + 5) 6 



(21) 



(22) 



which completes our calculation of ir n for all n. Figure 1 
shows the resulting QSD for ft = 100 and 7 = 1/2. 

Decay time. Having found the QSD we can calculate 
the decay time r. We return to Eq. (O and sum it up 
over n from zero to infinity. By virtue of Eq. (|13|) . the 
first term on the right tends to A s /(2fl), while the rest 
of the terms tend to zero. We obtain 



E : 



As_ 

2ft 



«*(! - S) 
2tt(1 + <5) 



-OAs 



(23) 
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FIG. 2: (color online). The quasi-stationary distribution 
(QSD) for fi = 100 and 7 = 1/2. The QSD includes four 
overlapping asymptotes: the fast-mode (1), the slow-mode 
(2), the boundary-layer (3), and the Poisson distribution (|18[) 
(triangles). The small mismatch between the curves is due 
to higher-order corrections. We found numerically that the 
mismatch goes down as ~ Sl -1 ^ 2 at large SI. 



which is exponentially small as long as HAs ^> 1. The 
decay time, in physical units, is therefore 



1 



2AI + S) gOAs 



1*6(1-8) 

As is monotone increasing with 5; its asymptotes are 



(24) 



As 



(2/3)£ 3 + (4/15 )<5 5 

2-Wy/2(l-5) + . 



., <5« 1, 



The exponent e nAs in Eq. (j24|) coincides with that ob- 
tained by Elgart and Kamenev Q who only considered 
the leading order of (a different version of) WKB. Our 
result (pM)) goes beyond the leading order and includes 
a pre-exponcnt. The pre-exponent diverges as S — > 1 
(when the immigration is stopped), so the decay time r 
diverges. This result could not have been predicted in 
the leading order of WKB 0. 

The assumptions made in the process of derivation of 
our results include the strong inequalities SI As » 1 , 
I2 << ni—n\ and l\ << min (n\, n%— n\). For 5 = 0(1), 
all the assumptions holds for f2 3> 1. For < 6 -C 1 
[above but close to the bifurcation point of the birth of 
the fixed points n± and rii of Eq. ([T])] the criterion is more 
stringent: S!(5 3 >• 1. 

In conclusion, by using a simple birth-death process 
as an example, we have developed a systematic the- 
ory of unlimited population growth driven by demo- 
graphic noise. We have found the complete metastable 
probability distribution of the population at long times, 
P n (t) ~ 7r„ e~*/ T , and the previously unknown important 
pre-exponential factor in the decay time r ~ e~ nAs . As 
7r„ has a power-law tail ~ n~ 2 , no distribution moments, 
except the zeroth one, exist. As a result, the initial- value 
problem for the master equation @ is highly singular. 



When starting from a well-behaved P n (t = 0), all of the 
distribution moments, except the zeroth one, diverge al- 
ready at t > 0. 

A general outcome of this work is that slow WKB 
modes should play an important role, along with fast 
WKB modes, in population escape problems. Consider, 
as an example, the Schlogl model [121 ] where, in addition 
to our three reactions, one also has 3A — ► 2A. For very 
small v the rate equation has an additional attracting 
point 77,3 3> «2- Here the stochastic population switches 
randomly between two metastable states peaked at n\ 
and 713. Now, if P n (t = 0) is located around n = m, 
there is an exponentially long intermediate regime when 
the probability flux is directed from n = n\ to n = 713, 
whereas the reverse flux is negligible. This important 
regime is accurately captured by the solution presented 
above. The slow-mode component of the solution (which 
describes deterministic motion "down the hill") is vital 
in determining the pre-exponents of the QSD and of the 
decay time. 

We acknowledge a useful discussion with Alex 
Kamenev. This work was supported by the Israel Sci- 
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